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SELF-ORGANIZED CRITICALITY, MULTIFRACTALS, AND 


INTERMITTENT TURBULENCE IN EARTH’S MAGNETOTAIL 
I. SYNOPSIS. 

We have achieved all the goals stated in our grant proposal in collaboration with Dr. 
C.C. Wu of the University of California at Los Angeles. Specifically, these include: 

• The understanding of the complexity induced nonlinear spatiotemporal structures and 
the coexisting propagating modes. 

• The development of plasma resonances and coherent structures in space plasmas. 

• The study of fluctuation-induced nonlinear instabilities and their role in the 
reconfiguration of magnetic topologies in the magnetotail. 

• The development of “anisotropic three-dimensional complexity” in the plasma sheet 
due to localized merging and interactions of the magnetic coherent structures and 
associated topological phase transitions. 

• The understanding of the intermittent turbulence and energization process of the 
observed Bursty Bulk Flows (BBF’s) in the Earth’s magnetotail. 

• The acceleration of ions due to the intermittent turbulence of propagating arid 
nonpropagating fluctuations 

In the following, we include lists of our published papers, invited talks, and 
professional activities. A detailed description of our accomplished research results is 
given in Section IV. 
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T. Chang, "Complexity" Induced Plasma Turbulence in Coronal Holes and the Solar 
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Transitions — The Origin of the Low Frequency Broad-Band Spectrum in the 
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2001. 
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III. ACTIVITIES. 


• Guest Co-Editor: Special Topics Section of Physics of Plasmas on “Modes of 
Transport", 2002. 

• Co-Convener: Special sessions on “Transport Phenomena in Space Plasmas”, AGU 
Fall Meeting in San Francisco, December 2001 . 

• Co-Convener: International Workshop on "Complexity in Earth's Magnetotail", 
Venice, Italy, 2002. 

• Co-Guest Editor: Special Issue on “Complexity in Earth's Magnetotail”, Annales 
Geophysicae, to be published in 2003. 

• Co-Convener, Symposium on Solar System Physics, Joint EGS-AGU-EUG 
Assembly, Nice, France, 2003. 

• Scientific Advisor: Goddard Space Flight Center [Group on Self-Organized 

Criticality]. 

• Member: Program Committee, International Workshop on Nonlinear Waves and 
Chaos in Space Plasmas, India, 2002. 

Invited Lectures: 

• Conference on Sun Earth Connection - Multiscale Coupling of Sim Earth System, Big 
Island, Hawaii, 2004. 

• Triennial Tutorial Lecturer for the 2002 General Assembly of the International Union 
of Radio Science, Maastricht, Netherlands. 

• Symposium on Complexity in Space Plasma Processes, IUGG Assembly, Sapporo, 
Japan, 2003. 


5 


• Per Bak Memorial Session on Nonlinear Processes in Space Plasmas, Joint EGS- 
AGU-EUG Assembly, Nice, France, 2003. 

• Session on Solar System Physics, Joint EGS-AGU-EUG Assembly, Nice, France, 
2003. 

• PARS (Polar Aeronomy and Radio Science) Workshop, Lake Arrowhead, CA, 2002. 

• Huuntsville Workshop on Particle Acceleration in Space and Astrophysical Plasmas, 
Huntsville, TN, 2002. 

• Eringen Symposium, Society of Engineering Science, State College, PA, 2002. 

• International Workshop on Complexity in the Magnetotail Dynamics, Venice, Italy, 

2002. 

Visiting Scientists and Students: 

Sandra Chapman [University of Warwick]. 

Nick Watkins [British Antarctic Survey]. 

George Rowlands [University of Warwick] 

Jay Johnson [Princeton Plasma Physics Laboratory]. 

Vincenzo Vitelli [Imperial College]. 

Andrew Yau [University of Calgary]. 

Tom March [University of Warwick]. 

James Merrifield [University of Warwick]. 

Listing, 2000 Outstanding Scientists of the 20th Century. 
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IV. Complexity induced anisotropic bimodal intermittent turbulence in space 
plasmas 

Tom Chang and Sunny W.Y. Tam 

Center for Space Research, Massachusetts Institute of Technology, Cambridge, 
Massachusetts 02139 USA 

Cheng-chin Wu 

Department of Physics and Astronomy, University of California, Los Angeles, California 
90095 USA 

Abstract 

The "physics of complexity" in space plasmas is the central theme of this exposition. 
It is demonstrated that the sporadic and localized interactions of magnetic coherent 
structures arising from the plasma resonances can be the source for the coexistence of 
nonpropagating spatiotemporal fluctuations and propagating modes. Non-Gaussian 
probability distribution functions of the intermittent fluctuations from direct numerical 
simulations are obtained and discussed. Power spectra and local intermittency measures 
using the wavelet analyses are presented to display the spottiness of the small-scale 
turbulent fluctuations and the non-uniformity of coarse-grained dissipation that can lead 
to magnetic topological reconfigurations. The technique of the dynamic renormalization 
group is applied to the study of the scaling properties of such type of multiscale 
fluctuations. Charged particle interactions with both the propagating and nonpropagating 
portions of the intermittent turbulence are also described. 
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I. INTRODUCTION 


Simultaneous coexistence of propagating modes and intermittent nonlinear 
spatiotemporal structures is the norm of the state of the plasma medium in the space 
environment. The "physics" of the bimodal state of such type of admixture of turbulent 
fluctuations may be understood from the point of view of the development and 
interactions of coherent structures arising from plasma resonance sites . 1 ' 6 

In this treatise, we shall consider the dynamical complexity in space plasmas from 
such a concept. Results of two-dimensional direct numerical simulations 7 ' 9 including the 
calculated fluctuation probability distribution functions and local intermittency measures 
based on the wavelet transforms will be presented to characterize the sporadic, localized, 
and scale-dependent nature of the intermittent turbulence. 

The concepts of scale invariance and symmetry-breaking phenomena of such 
invariance properties will be described based on the dynamical theory of the 
renormalization group . 10 Illustrative examples will be provided to elucidate the utility of 
this powerful theoretical technique in addressing the complexity of space plasmas. 

The intricate interactions among the charged particles in the plasma medium and the 
propagating and nonpropagating intermittent fluctuations will be considered and an 
example related to the energization of the auroral ions will be provided . 11 

We have endeavored to make this exposition sufficiently self-contained. 
Approximately one third of this treatise is a review of previous work and the remaining 
two thirds are discussions of new research results. 


8 



II. THE ORIGIN OF COMPLEXITY 


A. Alfvenic resonances and coherent structures 

Most field theoretical discussions begin with the concept of propagation of waves. For 
example, in the MHD formulation, one can combine the basic equations and express them in the 
following propagation forms: 

pcNldt = BVB+ , dB/dt = B'V\ + ■•• (1) 

where the ellipses represent the effects of the anisotropic pressure tensor, the compressible and 
dissipative effects, and all notations are standard. Equations (1) admit the well-known Alfven 
waves. For such waves to propagate the propagation vector k must contain a field-aligned 
component, i.e., B V— »/k-B*0. However, at sites where the parallel component of the 
propagation vector vanishes (i.e., at the resonance sites), the fluctuations are localized. Around 
these resonance sites (usually in the form of curves), it may be shown that the fluctuations are 
held back by the background magnetic field, forming Alfvenic coherent structures. 1-4 ’ 12 ’ 13 

B. Coarse-grained helicity 

Let us now consider the geometiy of the Alfvenic coherent structures. For an ideal MHD 
system, it has been suggested by Taylor 14 that in a relaxed state such a structure would be 
approximately force-free (i.e., JxB = 0) due to the approximate conservation of the coarse- 
grained helicity defined as K = JAB dV integrated over the coherent structure, where J and B 
are the current density and magnetic field and A is the vector potential. 

To obtain some physical insight of these structures, let us consider the special situation for the 
auroral region and/or the solar wind and make the reasonable assumption that the perturbed 
magnetic field fluctuations are much smaller than and essentially transverse to the mean magnetic 
field Bq (which will be temporarily assumed to be uniform for the current discussion). Thus, let 
us write B = ( SB x ,SBy,B 0 ) , where z is in the direction of the mean magnetic field, and (x, y) are 
orthogonal coordinates normal to z. The force-free condition for constant B 0 and V • J = 0 then 
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leads approximately to the scalar condition B- VJ Z = 0 , obtained by taking the z-component of 

the curl of JxB = 0 . n ’ 15 ’ 16 it can be shown that, with the inclusion of the kinetic effects through 
the anisotropic pressure terms and the generalized Ohm’s law, the above results are still 
approximately valid. We have, then, approximately, 

B 0 dJ z / dz = ~(SB x d / dx+SByd ldy)J z + • • • (2) 

where the ellipsis represents the other modifying effects. For convenience, let us introduce the 
flux function y/ by writing {d yrid y,-dy/id x) - (SB x ,SB y ) for the perturbed transverse 

components of the magnetic field in the (x,y) directions such that V B = 0 is satisfied. Then, 

J z and iff are governed by Eq. (2) and the Ampere’s law (neglecting the modifying effects 
represented by the ellipsis). 

A simple example of the flux function and axial current density satisfying the above 
conditions would be the class of circularly cylindrical solutions of if/(r) and J z (r) . Generally, 

the solutions would be more involved because of the variabilities of the local conditions of the 
plasma and the three-dimensional geometry. Moreover, the dynamic coherent structures with the 
inclusion of plasma pressure and other modifying effects (including electron-inertia terms) would 
be even more complicated. However, we expect these structures to be usually in the form of 
field-aligned flux tubes. Fig. 1 . 

Generally, there exist various types of propagation modes (whistler modes, electromagnetic 
ion cyclotron waves, etc.) in magnetized plasmas. Thus, we envision a corresponding number of 
different types of plasma resonances and associated coherent structures that typically characterize 
the dynamics of the plasma medium under the influence of a background magnetic field. 

Generally, such coherent structures may take on the shapes of convective forms, nonlinear 
solitary structures, pseudo-equilibrium configurations, as well as other types of spatiotemporal 
varieties. These structures might be locally generated or convected from elsewhere (such as some 
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of the observed structures in the solar wind that might have been originated from the Sun's 
surface 173 . Some of these structures may be more stable than the others. They, however, 
generally are not purely laminar entities as they are composed of bundled fluctuations of all 
frequencies. Because of the nature of the physics of complexity, it will be futile to attempt to 
evaluate and/or study the details and stabilities of each of these infinite varieties of structures; 
although some basic understanding of each type of these structures will generally be helpful in 
the comprehension of the full complexity of the underlying nonlinear plasma dynamics. 

These coherent structures will wiggle, migrate, deform and undergo different types of motions 
and interactions under the influence of the local plasma and magnetic topologies. In the next 
section, we will consider how the coherent structures can interact and produce the type of 
intermittency generally observed in a complex dynamical plasma. 

C. Interactions and complexity 

When coherent magnetic flux tubes of the same polarity migrate toward each other, strong 
local magnetic shears are created. Fig. 2. It has been demonstrated by Wu and Chang 7 ' 9 that 
existing sporadic nonpropagating fluctuations will generally migrate toward the strong local shear 
region. Eventually the mean local energies of the coherent structures will be dissipated into these 
concentrated fluctuations in the coarse-grained sense and induce reconfigurations of the magnetic 
field geometry. Figure 3(a) displays 2D MHD simulation results for homogeneous turbulence. 
The calculations were carried out with (512x512 ) grid points in a doubly periodic (x.y ) domain 

of length 2n in both directions. The left panel of Fig. 3(a) gives contours of the magnetic flux at 
t = 300 . For reference, the sound wave and Alfven wave traveling times through a distance 2 n 
are about 4.4 and 60, respectively. Regions of intense current density (strong magnetic shear) are 
evident in the right panel of Fig. 3(a). (We found the results are very similar with or without a z- 
component magnetic field.) Figure 3(b) is essentially the same as Fig. 3(a) except that these 
results were obtained with (1024x1024) grid points. Thus, our simulation results are quite 
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robust. To understand the effect of the coarse-grained dissipation, we have performed 2D MHD 
simulations of (512x512) grid points with a sheared magnetic field of B x = 0.1 cosy. Figure 

3(c) are the contours of the magnetic flux and current density distributions for the shear case. The 
neutral sheet regions are near y - Q.5n,\.5n. The clustering of the coherent structures (enhanced 
coarse-grained dissipation) and the reconfiguration of the contours of the magnetic fields near the 
neutral sheets are evident in the displays. 

Such enhanced intermittency at the intersection regions has been observed by Bruno et a/. 17,18 
in the solar wind using the tools of wavelet analyses and local intermittency measure (LIM). As 
mentioned above, the coarse-grained dissipation will then initiate "fluctuation-induced nonlinear 
instabilities"; 4,19 and, thereby reconfigure the topologies of the coherent structures of the same 
polarity into a combined lower local energetic state, eventually allowing the coherent structures to 
merge locally. On the other hand, when coherent structures of opposite polarities approach each 
other due to the forcing of the surrounding plasma, they might repel each other, scatter, or induce 
magnetically quiescent localized regions. Under any of the conditions of the above interaction 
scenarios, new fluctuations will be generated. And, these new fluctuations can provide new 
resonance sites; thereby nucleating new coherent structures of varied sizes. 

All such interactions can occur at any location of a flux tube along its field-aligned direction, 
and the phenomenon is fully three-dimensional. In order to gain some insight of the physical 
picture of the overall dynamics of the interactions of the coherent structures, let us again consider 
the auroral zone or the solar wind as an illustrative example. We make the plausible assumption 
that some aspects of the plasma dynamics may be approximately understood in terms of the 
formulation of reduced magnetohydrodynamics (RMHD). 20,21 In this approximation, we assume 
that the mean magnetic field is much larger than the transverse fields, and the field-aligned 
fluctuations of the magnetic and velocity components are much smaller than their transverse 
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counterparts. As a consequence, the density of the plasma is uniform. Writing the equations in 
SI units with p = 1 and M 0 = 1 , we have: 22,23 

dy/ / dt = B z d<p / dz , dco/dt-B-Vj (3) 

where yf(x,y) is the transverse flux function defined by B = e z xViff+ B z e z , <p(x,y) is the 
transverse stream function defined by v x = e z xV0, and co-^\<p is the vorticity, y = V 2 ^ is 
the field-aligned current density with d/dt-d/dt + \*V. The equations are written in the 
moving frame along the mean magnetic field direction z, and (x, y) are the transverse orthogonal 
directions. 

From Eqs. (3), we note that the primary nonlinear interactions occur generally in the 
transverse direction to the mean magnetic field. And, the coupling in the field-aligned direction 
is essentially linear. Thus, fluctuations generated by the transverse nonlinear interactions will 
scatter and evolve nonlinearly primarily in the transverse direction. At this point, we realize that 
the RMHD formulation is too restrictive as some of the interactions of the flux tubes may become 

more oblique and thereby allowing the fluctuations to attain a broader range of values of k\\ than 
otherwise would have been admitted by the RMHD approximation. Thus, a significant amount of 
the fluctuations generated by the interactions can become commensurate with the plasma 
dispersion relation and propagate in the field-aligned direction as Alfven waves due to this three- 
dimensional complexity-induced enhanced transport. Eventually a dynamic topology of a 
complex state of coexisting propagating and nonpropagating magnetic fluctuations is created. In 
the auroral region, the plasma may be electron-inertia dominated and the above discussion can be 
easily generalized to include such kinetic effects. 

HI. SCALING AND SYMMETRY BREAKING 
A. Path integral and the dynamic renormalization group 

In the above sections, we provided some convincing arguments as well as numerical and 
observational evidences indicating that space plasma turbulence is generally in a state of 
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topological complexity. By "complex" topological states we mean magnetic topologies that are 
not immediately deducible from the elemental (e.g., MHD and/or Vlasov) equations. 24 Below, 
we shall briefly address the salient features of the analogy between topological and equilibrium 
phase transitions. A thorough discussion of these ideas may be found in Chang 4,11,25 (and 
references contained therein). 

For nonlinear stochastic systems exhibiting complexity, the correlations among the 
fluctuations of the random dynamical fields are generally extremely long-ranged and there exist 
many correlation scales. The dynamics of such systems are notoriously difficult to handle either 
analytically or numerically. On the other hand, since the correlations are extremely long-ranged, 
it is reasonable to expect that the system will exhibit some sort of invariance under coarse- 
graining scale transformations. A powerful technique that utilizes this invariance property is the 
method of the dynamic renormalization group. 10,25,26 The technique is a generalization of the 
static renormalization group introduced by Wilson. 27 

As it has been demonstrated by Chang et al., 26 based on the path integral formalism, the 
behavior of a nonlinear stochastic system far from equilibrium may be described in terms of a 
"stochastic Lagrangian L ", such that the probability density functional P of the stochastic system 
is expressible as: 

-POpOM)) = J ex P {-/ jl( <p, 9, "l)didt } (4) 

where <p(x,/) = ^,(/ = 1,2,...,//) are the stochastic variables such as the fluctuating magnetic, 
velocity and electric fields, and X(x,t) = X, (/ = 1,2,..., V) are the conjugate stochastic momentum 
variables that may be rigorously derived from the underlying stochastic equations governing 
9 , 10,25,26 

Then, the renormalization-group (coarse-graining) transformation may be formally expressed 
as: 

dLld£ = RL (5) 
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where R is the renormalization-group (coarse-graining) transformation operator and £ is the 
coarse-graining parameter for the continuous group of transformations. It will be convenient to 
consider the state of the stochastic Lagrangian in terms of its parameters { P n } . Equation (5), 
then, specifies how the Lagrangian, L, flows (changes) with £ in the affine space spanned by 
{/>„}, Fig. 4. 

B. Forced and/or self-organized criticality 

Generally, there exists, a number of fixed points (singular points) in the flow field, at which 
dL! d£ = 0. At each such fixed point (L* or L** in Fig. 4), the correlation length should not be 
changing. However, the renormalization-group transformation requires that all length scales must 
change under the coarse-graining procedure. Therefore, to satisfy both requirements, the 
correlation length must be either infinite or zero. When it is at infinity, the dynamical system is 
then at a state of forced and/or self-organized criticality (FSOC), 25,28 analogous to the state of 
criticality in equilibrium phase transitions. 29 To study the stochastic behavior of a nonlinear 
dynamical system near such a dynamical critical state (e.g., the one characterized by the fixed 
point L*), we linearize the renormalization-group operator R about L*. The mathematical 
consequence of this approximation is that, close to dynamic criticality, certain linear 
combinations of the parameters that characterize the stochastic Lagrangian L will correlate with 
each other in the form of power laws. These include, in particular, the ( k , (o ), i.e. mode number 
and frequency, spectra of the correlations of the various fluctuations of the dynamic field 
variables. 

Such power law behavior has been detected in the probability distributions of solar flare 
intensities, 30 in the AE burst occurrences as a function of the AE burst strength, 31 in the global 
auroral UVI imagery of the statistics of size and energy dissipated by the magnetospheric 
system, 32 in the probability distributions of spatiotemporal magnetospheric disturbances as seen 
in the UVI images of the nighttime ionosphere, 33 and in the probability distributions of durations 
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of Bursty Bulk Flows ; 34 although some of the above interpretations of observed data may also be 
amenable to alternative explanations . 35 

In addition, it can be demonstrated from such a linearized analysis of the dynamic 
renormalization group that generally only a small number of (relevant) parameters are needed to 
characterize the stochastic state of the system near criticality 25 justifying the recent work 
suggesting that certain dynamic characteristics of the magnetotail could be modeled by the 
deterministic chaos of low-dimensional nonlinear systems . 36-38 
C. Illustrative examples 

The intermittency description for plasma turbulence of nonpropagating fluctuations may be 
modeled by the combination of a localized chaotic functional growth equation of a set of relevant 
order parameters and a functional transport equation of the control parameters. 6 ' 35 ' Below, we 
shall provide two simple phenomenological models, which may have some relevance, to the 
auroral zone or the solar wind . 5 
1. Model I 

Assuming that the parallel mean magnetic field Bo is sufficiently strong and the magnetic 
fluctuations dominate in the transverse directions, we introduce the flux function y/ for the 
transverse fluctuations as follows, 

B = e z xV V + B 0 e z (6) 

where we have set B z = B 0 = constant. The coherent structures for such a system are generally 

flux tubes approximately aligned in the mean parallel direction . 11 Conservation of helicity (e.g., 
under the RMHD approximation) indicates that the integral of y/ over a flux tube is 
approximately constant. Instead of invoking the RMHD formalism, however, here we simply 
consider yf as an order parameter. As the flux tubes merge and interact, they may correlate over 
long distances, which, in turn, will induce long relaxation times near FSOC . 5,25 Let us assume 
that the transverse size of the system is sufficiently broad compared to the cross sections of the 
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coherent structures (or flux tubes), such that we may invoke homogeneity and assume the 
dynamics to be independent of boundary effects. We may then model the dynamics of flux tube 
mergings and interactions, in the crudest approximation, in terms of the following order-disorder 
intermittency equation: 

dWk ldt=-Y k dFldW-k +fk (7) 

where V k are the Fourier components of the flux function, T* an analytic function of k 2 , 
F(y/k,k) the state function, and f k a random noise which includes all the other effects that are 

not included in the first two terms of this crude model. 

2. Modem 

In the above model, we have neglected both the effects of diffusion and convection. We next 
construct a phenomenological model that includes the transport of cross-field diffusion. We now 
assume the state function to depend on the flux function yf and the local pseudo-energy measure 
£ . Thus, in addition to the dynamic equation (7), we now also include a diffusion equation for 
£ . In Fourier space, we have 

d£ k /dt = -Dk 2 dF/d£- k +h k (8) 

where £ k are the Fourier components of £ , D(k) is the diffusion coefficient, and the state 
function is now F(yf k ,£ k ,k) , and h k is a random noise. By doing so, we separate the slow 
pseudo-energy transport due to diffusion of the local pseudo-energy measure £ from the noise 
term of (7). We note that an approach similar to these ideas have been considered by Klimas et 
a /. 40 

3. Dynamic renormalization-group analysis 

We have performed renormalization-group analyses as outlined above for the two kinetic 
models described above. We note that under the dynamic renormalization-group (DRG) 
transformation, the correlation function C of W k should scale as: 
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e a ^C\k,o)) = C(ke(,o)e aa/ ) (9) 

where co is the Fourier transform of t, £ the renormalization parameter as defined in Sec. IIIA, 
and (a c ,cko) the correlation and dynamic exponents. Thus, C/co a c /a<0 is an absolute invariant 

under the DRG, or C - ( o ~ * , where A = -a c / a m . DRG analysis of Model I with Gaussian noise 
yields the value of A to be approximately equal to 2.0. DRG analyses performed for Model II for 
Gaussian noises for several approximations yield the value for A to be approximately equal to 
1.88 to 1.66. 

Interestingly, for both models, DRG calculations give an approximate value of -1.0 for the 
co -exponent for the trace of the transverse magnetic correlation tensor. Matthaeus and 
Goldstein 41 had suggested that such an exponent might represent the superposition of discrete 
structures emerged from the solar convection zone. This value for the co -exponent is mean- 
field-like, thus is probably universally applicable to the low frequency fluctuations such as those 
considered by M & G as well as other more intermittent small scale fluctuations observed in 
space plasmas. Also, the corresponding ^-exponent is found to be approximately equal to -2 for 
both models. These results compare rather favorably with the results of our 2D MHD numerical 
simulations, Fig. 5. 

D. Symmetry breaking 

As the dynamical system evolves in time (autonomously or under external forcing), the state 
of the system (i.e., the values of the set of the parameters characterizing the stochastic 
Lagrangian, L) changes accordingly. A number of dynamical scenarios are possible. For 
example, the system may evolve from a critical state A (characterized by L**) to another critical 
state B (characterized by L *) as shown in Fig. 4. In this case, the system may evolve 
continuously from one critical state to another. On the other hand, the evolution from the critical 
state A to critical state C as shown in Fig. 4 would probably involve a dynamical instability 
characterized by a first-order-like topological phase transition (fluctuation-induced nonlinear 
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instability) because the dynamical path of evolution of the stochastic system would have to cross 
over a couple of topological (renormalization-group) separatrices. For such a situation the 
underlying magnetic topology and its related plasma state will generally undergo drastic changes. 
Similar ideas along these lines have been advanced by Sitnov et al n and simulated based on the 
cellular automata calculations of sandpile models. 43 ’ 44 

Under either of these above scenarios, the spectra indices will generally change either 
continuously or abruptly. Such type of multifractal phenomena is commonly observed in the 
magnetotail, the auroral zone, and the solar wind. 11,17 ’ 18,45 ' 49 

Alternatively, a dynamical system may evolve from a critical state A to a state D (as shown in 
Fig. 4) which may not be situated in a regime dominated by any of the fixed points; in such a 
case, the final state of the system will no longer exhibit any of the characteristic properties that 
are associated with dynamic criticality. As another possibility, the dynamical system may deviate 
only moderately from the domain of a critical state characterized by a particular fixed point such 
that the system may still display low-dimensional scaling laws, but the scaling laws may now be 
deduced from straightforward dimensional arguments. The system is then in a so-called mean- 
field state. (For general references of symmetry breaking and nonlinear crossover, see Chang and 
Stanley 50 ; Chang et al . 51,52 ; Nicoll et al . 53,S4 ) 

Experimental observations of plasma fluctuations in the space plasma environment generally 
yield broken power law spectra similar to those displayed in Fig. 5 of the 2D numerical 
simulation results. Such abrupt changes of scaling powers of the k-spectra are signatures of 
symmetry breaking. The broken symmetries may be due to the abrupt change of the degree of 
intermittency of fluctuations from large to small scales, or due to the change of the underlying 
physics (e.g., from MHD to kinetic processes), or variations of external forcing, or finite 
boundaries and other effects. 
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IV. DYNAMICAL INTERMITTENCY 


Nearly all fluctuations in space plasmas exhibit intermittency. For turbulent dynamical 
systems with intermittency, the transfer of energy (or other relevant scalars and tensors) due to 
fluctuations from one scale to another deviates significantly from uniformity. A technique of 
measuring the degree of intermittency is the study of the departure from Gaussianity the 
probability distribution functions of turbulent fluctuations at different scales. To demonstrate this 
point, let us refer to the 2D numerical simulation results described in Sec. II. For example, we 

may generate the probability distribution function P(SB 2 ,S) of SB 2 (x,S ) 

= B 2 (jc + S) - B 2 (x) at a given time t for such simulations, where S is the scale of separation in 
the x-direction. The left panel of Fig. 6(a) displays the calculated results of P(SB 2 ,S) from the 

(5 1 2x5 1 2) 2D simulation for the homogeneous case for several scales S . 

From this figure, we note that the deviation from Gaussianity becomes more and more 
pronounced at smaller and smaller scales. In an interesting paper by Hnat et al. 55 , they 
demonstrated that such probability distributions for solar wind fluctuations exhibit approximate 
mono-power scaling according to the following functional relation: 

P(8B 2 ,S) = S- s P s {SB 2 S- 5 ,S) (10) 

where s is the mono-scaling power. We found that mono-power scaling also holds approximately 
for our simulated results with the value of s equal to approximately 0.335 as shown in the right 
panel of Fig. 6(a). 

The reason for mono-power scaling for SB 2 may be understood in terms of the 

renormalization-group arguments presented in Sec. III. If we assume that SB 2 is one of the 
relevant eigenoperators near a critical fixed point, then the probability distribution function for 
P(SB 2 ,S ) , SB 2 , as well as 5 will scale linearly as follows: 

P' = P exp (a p l) , SB 2 ' = SB 2 exp(a /) , S' = 5 exp (a s l) (11) 
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where (a p ,a B i ,a s ) are scaling powers. Thus, we obtain two irreducible absolute invariants: 

p j /a s an( j fip 2 / b 2 since P — P{8B 2 ,8 ) , there must be a functional relation between 

these two invariants. 51,52 Therefore, we obtain the following scaling relation among 

(P,SB 2 ,S): PI8°' la ° =F{SB 2 / 8°*'°*). 

Without loss of generality, we may choose a s = 1. With the additional constraint that the 

probability distribution functions are normalized, we immediately obtain the expression of Hnat 
et al. as shown in Eq. (10). 

Actually the scaling relation (10) is approximate in that the tails of the scaled distributions in 
the right panel of Fig. 6(a) do not exactly fall onto one curve. This is the intrinsic nature of the 
strong intermittency at small scales. Thus, representations of the probability density functions 
will involve multi-parameters in general and may sometimes be represented more accurately, for 
example, by the Castaing or Kappa distributions 56 * 59 Their scaling properties are more subtle and 
will not be considered in this treatise. 

As it has been pointed out in Sec. IIC, sporadic nonpropagating fluctuations will generally 
migrate toward regions of strong'shear. We have calculated the probability distribution functions 
for regions near and away from the neutral sheet for the (512x512) 2D simulation for the shear 
case. Figs. 6(b,c). As expected, the degree of intermittency (deviation from Gaussianity) is much 
stronger near the neutral sheet than away from it. Both regions exhibit approximate mono-power 
scaling; though, again the tail regions do not scale exactly, particularly near the neutral sheet. 

Since the degree of intermittency generally increases inversely with scale, it will be interesting 
to study the degrees of intermittency locally at different scales. This can be accomplished by the 
method of Local Intermittency Measure (LIM) using the wavelet transforms. A wavelet 
transform generally is composed of modes which are square integrable localized functions that 
are capable of unfolding fluctuating fields into space and scale. 60 Left panel of Fig. 7(1) is the 


21 



power spectrum of a complex Morlet wavelet transform of the current density J z for the 
(5 1 2x5 1 2) 2D homogeneous case. We notice that the intensity of the current density is sporadic 
and varies nonuniformly with scale. 

We now define LIM(l) as the ratio of the squared wavelet amplitude |(/(x,£)| 2 and its space 

averaged value < ri>')| 2 > x . We note that LIM(l) = 1 for the Fourier spectrum. To 

emphasize the variation of intensity with scale, we also consider the logarithm of LIM(l). It has 
been suggested by Meneveau 61 that the space average of the square of LIM(l), which is a scale 
dependent measure of the kurtosis or flatness, is a convenient gauge of the deviation of 
intermittency from Gaussianity. We denote this measure by LIM(2). It is equal to 3 if the 
probability distribution is Gaussian. Bottom right panel of Fig. 7(1) and Fig. 8(a) are graphical 
displays of the calculated results of logLIM(l) and LIM(2) for the (512x512) 2D simulation for 
the homogeneous case using the complex Morlet transform. We notice that the fluctuations are 
indeed scale dependent, localized and strongly intermittent at small scales. Similar experimental 
results using the wavelet transforms have been found, for example, by Consolini et al} 2 for the 
magnetotail and Bruno et al} 1 for the solar wind. 

Figure 7(2) gives graphical displays of logLIM(l) for the (512x512) 2D simulation for the 
shear case for both away from and at the neutral sheet region. We note that the degree of 
intermittency is indeed much stronger at the neutral sheet region. This result is also confirmed by 
the graphical displays of LIM(2) in Fig. 8(b). 

In the above, we considered some simplified models and numerical examples that may have 
some relevance to intermittent turbulence in space plasmas. Realistic models for these 
phenomena will generally be much more complicated. For example, from the RMHD 
formulation of (3), we recognize that there should at least be two competing order parameters. 
These are the flux function \f/ and the stream function tp (which is linearly proportional to the 
electrostatic potential). Thus, the intermittency equation (such as (7)) needs to be generalized to 
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accommodate these coupled order parameters. Within the RMHD formulation, there exist useful 
Hamiltonian and operator-algebra structures, 63 which should prove invaluable in developing the 
generalized state function F of the coupled order parameters and the state variables as well as the 
intermittency equation itself. Formulations of coupled order parameters and their related 
theoretical analyses for a variety of criticality problems in condensed matter physics have been 
considered by Chang et al. }0 (and references contained therein). 

In addition, the transport equation such as Eq. (8) for the global system should generally also 
include convection and acceleration terms in addition to that of diffusion. 6 Thus, at the minimum 
our model transport equation must take on the form of the RMHD (3) with the addition of 
"coarse-grained" dissipation terms which generally will be functionals of the coupled order 
parameters. It should also contain terms representing the complexity-induced enhanced field- 
aligned transport. These generalizations will not be considered in this treatise. 

V. ENERGIZATION OF IONS BY INTERMITTENT FLUCTUATIONS IN THE 
AURORAL ZONE 

It has long been recognized that the commonly observed broadband, low frequency electric 
field fluctuations are responsible for the acceleration of oxygen ions in the auroral zone. In order 
for the fluctuating electric field to resonantly accelerate the ions continuously as the ions evolve 
upward along the field lines, they must be in continuous resonance with the ions. There did not 
seem to exist a fully viable mechanism that can generate a spectrum of fluctuations broadband 
and incoherent enough to fulfill this stringent requirement. 

Assuming that the RMHD formulation holds approximately in the auroral zone, the 
electrostatic fluctuations transverse to the field-aligned direction are given approximately by the 
velocity fluctuations: \xB z e z . The ordering due to the stream function <p may be important in 

the auroral zone and therefore, the electrostatic fluctuations can be quite significant there. 
Because of the small scales involved, the dynamic intermittency produced by the merging and 
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interactions of the coherent structures is probably generated by the whistler turbulence, electron- 
inertia related tearing modes, and/or other collisionless modes. 11 Therefore, a significant portion 
of the fluctuations would be kinetic. Nevertheless, the electric field fluctuations would still be 
predominantly transverse and electrostatic. Thus, a significant portion of the low frequency 
fluctuations commonly observed in the auroral zone are probably contributed by these 
nonpropagating intermittent fluctuations. 

Below, we shall briefly discuss how such fluctuations can efficiently energize the oxygen ions 
from ionospheric to magnetospheric energies. Assuming the oxygen ions are test particles, they 
would respond to the transverse electric field fluctuations E x near the oxygen gyrofrequency 
locally according to the Langevin equation: 

dv 1 /dt^q i E 1 /m i (12) 

To understand the stochastic nature of the Langevin equation, we visualize an ensemble of 
ions /(v ± ) and study its stochastic properties. Assuming that the interaction times among the 
particles and the local electric field fluctuations are small compared to the global evolution time, 
we may write within the interaction time scale: 

f(\ ± ,t + At) = -A\±,t)P,(\ 1 - Avi,Avj.)c/Avi (13) 

where P,(v 1 -Av±,Av 1 ) is the normalized transition probability of a particle whose velocity 

changes from v ± - Av ± to Vj. in At, and A\± ranges over all possible magnitudes and 
transverse directions. Standard procedure at this point is to expand both sides of (13) in Taylor 
series expansions: 

|^A, + O((40’) = -^((Av x )/] 

: [( Av i Ay i)/] + °(( Ay i)’) < 14 > 
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where 


( Av l)= J^(v x ,Av x )Av x rf(Av x ),and (15a) 

(Av x Av x ) = J/J(v 1 ,Av ± )Av ± Av ± rf(Av 1 ). (15b) 

If we assume the 0((Av x ) 3 ) terms are of order (At) 2 or higher, then in the limit of At — » 0 , we 
obtain a Fokker-Planck equation, 64,65 where the drift and diffusion coefficients are defined as: 
D] =< Avi > /At and D 2 =< AvjAvj. > /2At in the limit of At — » 0. These coefficients may 
be calculated straightforwardly using the Langevin equation. If the transition probability P, is 
symmetric in Av_i_ , then D| vanishes and (14) leads to a diffusion equation in the transverse 

direction. We note that if the electric field fluctuations are Gaussian, then the higher order 
correlations of the fluctuations are automatically equal to zero. 

We shall come back to the discussion of the effects of general intermittent fluctuations on 
particle energization processes. For the moment, let us assume the approach using the Fokker- 
Planck formulation is valid and proceed. Since we have assumed the time scale for the particle- 
fluctuation interactions is much smaller than the global evolution time of the ion populations, we 
may then write the steady global evolution equation along an auroral field line s under the 
guiding center approximation and neglecting the cross-field drift as: 66-68 

_a_r / 1 1 a r vi dB z \ 3 r ^ / ~ 

ds 11 B 2 2B Z ds B z v x dv x x 2B Z dsB z 


_L_a_ 

v x dv x 


n Z> ±1 


( 16 ) 


where V||,v x are the parallel and perpendicular components of the particle velocity with respect to 
the field-aligned direction. This expression may be interpreted as a convective-diffusion equation 
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for the density of the guiding center ions per unit length of flux tube f I B z , in the coordinate 
space of (s, V||, v_l). 

To evaluate D x , the gyrotropic perpendicular diffusion coefficient, from D 2 , we recognize 

that the electric field fluctuations are broadband both in k x and CO . Therefore, at all times, some 
portion of the fluctuations will be in resonance with the ions. The resonance condition, however, 
is strongly dependent on the localization and scale dependency of the intermittent fluctuations." 
We demonstrate below, as a simple illustrative example, how such resonant interactions may be 
accomplished by neglecting the Doppler shifts due to k such that only the intermittent 
fluctuations clustering around the instantaneous gyrofrequency of the ions provide the main 
contributions to the diffusion process. Standard arguments then lead to the following expression 
for the perpendicular diffusion coefficient: 

Z> 1 =(* ? , ! /2m/)(|£ ! |(Q,)) r (17) 

where < |£ 2 1 (£2,) > r is the resonant portion of the average of the square of the transverse electric 

field fluctuations evaluated at the instantaneous gyrofrequency of the ions, £2/ . 

Measurements by polar orbiting satellites indicate that the electric field spectral density Z 
follows an approximate power law lT a in the range of the local oxygen gyrofrequencies, where 
a is a constant. If we make the additional approximations by assuming that the spectrum 
observed at the satellite is applicable to all altitudes and choosing the geomagnetic field to scale 

with the altitude as s ~ 3 , we would then expect Z(£2, , 5 ) to vary with altitude s as s ia . Because 
we have made some rather restrictive resonance requirements for the fluctuations to interact with 
the ions, we expect the resonant portion of the average of the square of the transverse electric 
field fluctuations to be only a fraction TJ of the total measured electric field spectral density. 
Therefore, we arrive at the following approximate expression for the diffusion coefficient: 
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Z) i =( W , 2 /2 Wl 2 )Z 0 (5/5 o ) 3a (18) 

We have performed global Monte Carlo simulations for Eqs. (16) and (18) for the conic event 
discussed by Retterer et al. 61 with or = 1.7 and 2 0 = 1.9x1 O' 7 (V/m) 2 sec/ rad . Top panel of 
Figure 9 shows the measured oxygen velocity distribution contours and the second panel of Fig. 9 
shows the corresponding calculated contours for 7 ] = 1 / 8 at the satellite altitude of s 0 = 2 R E ■ 
Thus, with one eighth of the measured electric field spectral density contributing, the broadband 
electric field fluctuations can adequately generate an oxygen distribution function with the energy 
and shape comparable to that obtained from observations. We have also calculated the oxygen 
ion distributions for a range of altitudes under the same conditions. Figure 10 is a plot of the 
average parallel energy versus the average perpendicular energy per oxygen ion as the ions 
evolve upward along the geomagnetic field line. We note that as the energies increase with 
altitudes, the ratio of the energies becomes nearly a constant. 

These results are comparable to our previous calculated results based on the assumption that 
the relevant fluctuations were purely field-aligned propagating electromagnetic ion cyclotron 
waves. 66-68 As discussed in the previous sections, we generally expect the coexistence of 
nonpropagating transverse electrostatic nonlinear fluctuations and a small fraction of field-aligned 
propagating waves in the auroral zone. Thus, the ion energization process in the auroral zone is 
probably due to a combination of both types of fluctuations. As it has been discussed in Chang et 
al. 66 , an asymptotic solution exhibiting such behavior may be obtained analytically in closed 
form. Therefore, in the asymptotic limit (i.e, at sufficiently high altitudes), it is expected that 
such an ion distribution will become entirely independent of its low altitude initial conditions. In 
fact, it has been shown by Crew and Chang 68 that the ion distributions will become self-similar at 
sufficiently high altitudes and eveiything will scale with the altitude. 

The above sample calculations did not include the self-consistent electric field that must be 
determined in conjunction with the energization of the ions as well as the electrons. This is 
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particularly relevant in the downward auroral current region where the electric field can provide a 
significant pressure cooker effect such as that suggested by Gomey et al. 69 and demonstrated 
convincingly by Jasperse 70 and Jasperse and Grossbard 71 based on global evolutional calculations 

•t'l nr nf a *1*1 

similar to those considered by Tam and Chang * for the solar wind and Tam et al. ’ for the 
polar wind. These ideas will not be considered in this treatise. 

We now return to the discussion of the effect of intermittency on ion heating. Measurements 
of the electric field spectral density are generally limited by the response capabilities of the 
measuring instruments. The faster the instruments can collect data, the more refined the scales of 
the measurements. As it has been seen in Sec. IV, we expect the measured spectrum density to 
exhibit small-scale intermittency behavior. In fact, it is known that fast response measurements 
generally exhibit strongly intermittent signatures of the fluctuations. In the diffusion 
approximation, the ion energization process is limited by the amplitude of the second moment of 
the probability distribution of the fluctuations. This amplitude may become smaller as the scale 
of measurements is reduced. Thus, in the limit of small scales, the amplitude of the measured 
spectrum may decrease and thereby requiring a larger value of rj to accomplish the same level of 
energization. 

But, the effects of the intermittency of the fluctuations on particle energization may be 
underestimated if we stay within the diffusion approximation. As it can be seen from the 
derivation of the diffusion approximation above, only the second order correlations of the 
fluctuations were included in the energization process. Since for intermittent turbulence, the 
probability distributions of the fluctuations are generally non-Gaussian, the effects of the 
intermittency can manifest in the higher order correlations beyond the second order diffusion 
coefficient. This implies that the higher order correlations of the velocity fluctuations may be of 
the order of At and therefore cannot be neglected in Eq. (14). Under such circumstances, the 
Fokker-Planck and diffusion approximations of the ion energization processes can become 
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inadequate. A more appropriate approach to address such non-Gaussian stochastic processes is to 
refer directly to the functional equation (13) using the non-Gaussian transition probability or the 
Langevin equation (12) with the actual intermittent time series of the electric field fluctuations. 




We have performed global simulations based on Eq. (13) for non-Gaussian intermittent 
fluctuations exhibiting the shape suggested by Castaing et a/. 58 : 

In 2 (cr/p-p)^ 


H^ ) = il eXP 


/ tl > 

G 

f 


exp 

2cr J 

V 


2A 2 


da 

J 


(19) 


where £ represents either the x- or ^-component of the dimensionless transverse velocity 


fluctuations and A> 0 is a parameter that characterizes the intermittency. We set In a$ = -A 2 , 

to ensure the variance equal to unity. For A- 0 , Eq. (19) reduces to a Gaussian distribution. As 
A increases, the degree of intermittency increases. Third panel of Fig. 9 shows the contours 
calculated for A = 1 with 1) = 1 / 8 . For this case, the degree of intermittency is not strong 
enough to significantly affect the value of 77 . But with strong intermittent fluctuations (A = 2, 
bottom panel of Fig. 9), a value of tj equal to 1/5 is required to adequately generate the ion conic 
to observed energies. 

VL SUMMARY 

We have provided a modem description of dynamical complexity relevant to the intermittent 
turbulence of coexisting nonpropagating spatiotemporal fluctuations and propagating modes in 
space plasmas. The theory is based on the physical concepts of sporadic and localized 
interactions of coherent structures that emerge naturally from plasma resonances. 

The technique of the dynamic renormalization group is applied to the study of forced and/or 
self-organized criticality (FSOC), scale invariance, and symmetry breaking, related to such type 
of multiscale fluctuations. We also demonstrated that the particle interactions with the 
intermittent turbulence could lead to the efficient energization of the plasma populations such as 
auroral ions. Numerical examples are presented to illustrate the concepts and methodology. 
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Figure Captions 

1 . Field-aligned spatiotemporal coherent structures. 

2. Cross-sectional view of coherent structures of the same polarity. Contours indicate constants 
of magnetic flux and arrows indicate directions of magnetic field. Blackened area is an 
intense current sheet (strong magnetic shear). 

3. 2D MHD simulations of interacting coherent structures. Left panels are contours of 
magnetic flux at / = 300 . Corresponding current density distributions are given in the right 
panels. For reference, the sound wave and Alfven wave traveling times through a distance 
2 n are about 4.4 and 60, respectively. 

4. Renormalization-group trajectories and fixed points. 

5. Fourier spectra of B 2 (k) at / = 300 (solid curves) and 600 (dashed curves) for' the 2D 
MHD simulations with (512x512) and (1024x1024) grid points. The straight lines have 
slopes of -2. Note that the scaling range is more extended for the larger simulation. 

6. 2D MHD simulation results with (512x512) grid points for the probability density 
distributions of the difference of the square of the magnetic field fluctuations at scales of S = 
2 (red), 4 (green), 8 (black), 16 (blue) and 32 (magenta) units of grid spacing. Left panels 
display the calculated results of P(SB 2 ,S ) and the right panels are the corresponding 
scaled distributions. 

7. Wavelet analyses of the 2D MHD simulation results with (512x512) grid points. Top 
panels show raw spatial distributions of J z or B . Bottom left panel of (1) gives the power 
spectrum of complex Morlet wavelet transform of the current density J z and right bottom 

panel of (1) gives the contours of the corresponding logLIM(l) of B for the homogeneous 
case. Left bottom panel of (2) gives the contours of logLIM(l) of B 2 at the neutral sheet 
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region and right bottom panel gives the corresponding contours away from the neutral sheet 
region for the shear case. The x-axis and scales are in units of the grid points. 

8. L1M(2) of B 2 using the Morlet wavelet transform for (a) the homogeneous case and (b) 
shear case (+++ for y = 1 .5/r and °°° for y = 1 ). 

9. Observed and calculated velocity contour plots for conic event of Retterer et al. 61 Top panel: 
observed contours. Second panel: Simulation results based on the diffusion approximation 
( A = 0). Third panel: Simulation results for weak intermittency (A = 1). Bottom panel: 
Simulation results for strong intermittency (A = 2). 

10. Solid line depicts W\\ versus W ± for simulated conic events. Dashed line is the asymptote 
predicted by Chang et al ^ 
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Fig. 1 
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(a) Homogeneous Case (512 x 512) 


(b) Homogeneous Case (1024 x 1024) 


(c) Shear Case (512 x 512) 
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Fig. 7 

(1) Homogeneous Case 


a) J z along y=n a) B 2 along y= 71 



(2) Shear Case 


a) B 2 along y= 1 .5n a) B 2 along y= 1 .Bn 
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